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Abstract 

A numerical approach to disordered 2D superconductors described by BCS 
mean field theory is outlined. The energy gap and the superfluid den- 
sity at zero temperature and the quasiparticle density of states are studied. 
The method involves approximate self-consistent solutions of the Bogolubov- 
de Gennes equations on finite square lattices. Where comparison is possible, 
the results of standard analytic approaches to this problem are reproduced. 
Detailed modeling of impurity effects is practical using this approach. The 
range of the impurity potential is shown to be of quantitative importance in 
the case of strong potential scatterers. We discuss the implications for exper- 
iments, such as the rapid suppression of superconductivity by Zn doping in 

Copper-Oxide superconductors. 
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I. INTRODUCTION 



It is well established that chemical substitutions on copper oxide superconductors can 
qualitatively alter the properties of these materials in both their normal and superconducting 
states^!. In the normal state,impurity substitution can be used to test models of electronic 
transport. In the superconducting state, impurity scattering effects are sensitive to order 
parameter symmetryi as well as other properties. It is well known that non-magnetic impu- 
rities are pair breakers in d-wave or other non-trivial pairing states with nodes of the energy 
gap on the Fermi surface. They produce a finite lifetime of the quasiparticles around the 
gap nodes and a finite density of states at low energy. 

In order to parameterize impurity scattering, a physical picture for the normal state, 
the superconducting state and the impurity is required. In this paper we are motivated 
by Zn doping experiments in cuprates which we assume can be treated as a pure potential 
scatterer of fermions which obey Luttinger's theorem. Since d-orbitals of Zn ++ are fully 
occupied, it is naively expected to behave as a non-magnetic impurity. However, electronic 
correlations may modify this picture somewhat. Nuclear Magnetic Resonance measurements 
by Mahajan and co-workersi show that the copper spin correlations are severely modified on 
neighboring Cu sites, which could in principle give rise to effects analogous to the spin-flip 
scattering process of conventional magnetic impurities. As discussed recently by Borkowski 
and HirschfeldS, the unknown relative strength of spin flip and impurity scattering rates in 
Zn doped systems poses an obstacle to the quantitative analysis of experiments involving 
Zn even within the conventional BCS formalism. 

Even if spin-flip scattering is negligible, the modification of copper spin correlations 
in the vicinity of an impurity site has the effect of increasing the range of an effective 
scattering potential. Information about the effective potential in Y Ba%{Gu\- x Zn^)iP%M 
can be obtained from the residual sheet resistance estimated by extrapolation from the 
normal state. The residual resistance from two dimensional potential scattering is given by 
the approximate expression! 
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Si are the phase shifts which are constrained to satisfy the Friedel sum- rule AZ = - 
AZ is the difference in the number of conduction electrons in the system without impurities 
and with one impurity. If we assume Zn removes one electron from the conduction band 
(i.e. AZ = — 1) and the effective potential is an impenetrable disc of radius a, we estimate 
that the residual sheet resistance is ~ O.QkQ for x = 0.01 in the I = dominant scattering 
channel (a ~ 0.54/ kf), which is a factor three smaller than the experimental value 
This discrepancy may indicate appreciable phase shifts in one or more higher angular mo- 
mentum channels. In fact the I = 2 "near resonant channel" , which requires a larger radius 
of the scattering potential (a ~ 2.8/kfp yields ~ 1.8M1 for x = 0.01. Recently Poilblanc, 
Scalapino, and Hankei have investigated the effects of non-magnetic impurities in antiferro- 
magnetically correlated systems. They find that the scatterings in I < 2 channels are strong 
in that system. 

The above discussion suggests that detailed structure of the impurity potential is impor- 
tant in making a quantitative study of Zn doped materials. As seen below, lattice effects may 
also be of quantitative importance for short coherence length superconductors. Recently, 
in a brief reporti, we presented a numerical study of the disorder effect in two dimensional 
superconductors of various pairing symmetries in the limit of strong impurity potentials. 
We demonstrated that a short but finite range potential has a much stronger effect than 
an on-site ( "5- function" ) potential. In fact the finite range of the potential may be the pri- 
mary reason for the rapid suppression of pairing correlations with impurity concentration in 
YBa(Cu, Zn)0. The importance of finite potential range has been pointed out recently by 
Balatsky et al0 in connection with non-universal behavior of the low frequency conductance 
in d-wave superconductors. 

In this paper we present a more detailed study of the disorder effect in two dimensional 
superconductors. The model used here is a lattice BCS mean-field Hamiltonian with disorder 
defined by, 



H[A rT ] = -t 4*Cr'o + £(Arr4 T 4+r| + h.C.) + £(£ r - //)4 CT C rCT , (2) 

(rr')cr rr rcr r, 

where (rr') denotes nearest neighbors, \i is the chemical potential. V Yi)T is a scattering 
potential of an impurity at rj. Lacking detailed knowledge of the scattering potential in 
high-T c cuprates, we assume a model potential form for it: V rur =Vo5 ri!r + Vi(5 ri - r ±x + 
5 ri ^ r ±y). When V\ = 0, it is a 5-function potential, otherwise it is finite ranged. A rT is the 
superconducting gap parameter. For a given filling factor of electrons n e , A rr and \i should 
be determined self-consistently from the relations 

A rr = J(c rT c r+Ti ) ArT , (3) 

U e = ^£(4crCrcr)A rT , (4) 
rcr 

where J is the coupling constant (assumed disorder independent), N is the system size, and 
{A)A rT means the average of A in the ground state of H[A rT \. Without disorder A rr is 
independent of r and has a particular symmetry with respect to r. In this paper only the 
on-site s-wave pairing state A rr = AS T o and the d-wave pairing state A rT = A(5 T ± X — S T ±y) 
in two dimensions will be considered. 

The Hamiltonian (|2]) is bilinear in fermion operators and can be diagonalized by solving 
a one-particle problem. A rr and ji should be determined self-consistently from Eqs. |3| 
and |], which amount to solving the Bogolubov-de Gennes equations for the disordered 
superconductor. For simplicity in calculation, we perform a particle-hole transformation for 
the down-spin electrons, i.e. c r | « — > cL, and re-express (|2|) as 

H'[A TT ] = -t £ (4 T C r /| - cJlCru) + £( A rrC rT C r+ri + h.C.) 
(rr'> rr 

+ Z)E V n,r ~ V) (4t c rt - 4| c r|) + const. (5) 

r r, 

H' has the usual tight binding model form, but the hopping constant, the chemical potential, 
and the impurity potential have opposite signs for the up spin electrons and the down spin 
electrons0. In the remainder of the paper we set the hopping constant t=l. 

In Sec. II, the numerical method and essential approximation is outlined. Results for 



the disorder dependence of the zero temperature gap and superfTuid density p s as well as 
the quasiparticle density of states p are presented. In Sec. Ill, a concluding remark is given. 

II. NUMERICAL METHOD AND RESULTS 
A. Gap parameter and superfluid density 

In the presence of disorder, the energy gap is space dependent. We determine A rr by it- 
eratively solving H' with the self-consistent conditions. We start from an initial gap function 
A rr with a certain pairing symmetry. After diagonalizing H', we find a new gap function A rr 
from the self-consistent equations and then use it as input to repeat the above process until 
the self-consistent conditions are satisfied. This is a strict self-consistent iterative process. 
However, since the gap function at every site needs be adjusted to satisfy the self-consistent 
equations, this is excessively time consuming when an average over a large number of impu- 
rity configurations is required. We shall however perform this strict self-consistent iteration 
procedure only for studying the properties of a single impurity system and for checking the 
accuracy of the approximation used in many-impurity cases. 

As mentioned above, to solve the self-consistent equations, the Hamiltonian needs be 
exactly diagonalized. This can be done, however, only on small lattices in the presence of 
disorder. For a superconductor, a characteristic length scale is the superconducting corre- 
lation length £ ~ hvp/irA. If £ is larger than the dimension of the system, the finite size 
effect is large and the analysis of the disorder effect may be subtle. To avoid this situation, 
we shall limit our calculations only to cases where £ is much smaller than the dimension of 
the system. 

The existence of a finite superfluid density p s is a defining property of superconductors. 
Experimentally p s is determined from the microwave measurement of the penetration depth. 
p s on a finite lattice can be evaluated directly from the current-current correlation function, 
using the eigenfunctions obtained from the iteration procedure described above,0 
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where (A) = Tr(Ae~ pH ) /Tre' 1311 ', (i^) is the kinetic energy along x-direction, and 



Before studying disordered systems, we consider a one-impurity system. We first consider 
the change of the gap function induced by an impurity. Fig. p] shows the self-consistent 
energy gap A rT for a s-wave superconducting state with one impurity on a 21x21 lattice. 
Since the scattering potential is short ranged, the gap function changes only in the vicinity 
of the impurity. Far away from the impurity, A rT approaches to the value of the energy gap 
without disorder. In a region with a length scale of the range of scattering potential around 
the impurity site A r is largely reduced due to the strong suppression of the probability of 
an electron hopping to this region by the impurity scattering. Beyond this region, there 
exists a relatively larger region with a length scale comparable with the superconducting 
correlation length £ where a weak oscillation of A rT in space is observed. This oscillation is 
due to the interplay between the impurity scattering and the superconducting correlations. 
Because of the lattice effect, A rr is not isotropic in space. It is relatively strong along the 
diagonal direction. Along other directions, it is too weak to be resolved from the figure. For 
the d-wave state, similar results have been found. But the oscillation of A rT along two axes 
seems more apparent in this case. 

For finite impurity doping systems where many impurity configurations are used in the 
disorder average we shall approximate A rT in (0) by the average of the pairing correlation 
functions (c r jC r+T |) in space, A T , so that only a simplified self-consistent equation, 



with c„ = Er CT 0n( 







A r = TfH ^(Cr|C r+rT )A r 



(9) 



r 
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needs be solved. This approximation is to ignore the fluctuation of A rT (but not (c r jC r+r |)) 
in space. 

The errors resulting from the above approximation can be found by directly comparing 
the results obtained with and without the approximation. We have calculated the errors 
for the pairing amplitude for several arbitrarily chosen configurations of impurities in both 
strong and weak scattering potential limits. For all the cases we have studied, we find that 
the relative errors in the average energy gap are small compared with the combined errors 
produced by the disorder average or the finite size effect. For example, for a randomly chosen 
system of 3 impurities on a 14x14 lattice, the relative errors in the average energy gap are 
0.1% (0.6%) in a weak potential V = 2 and V\ = and 2%(2%) in a strong potential V = 20 
and V\ = for the s-wave (d-wave) pairing state. As the impurity concentration increases, 
the error increases slightly. The errors for the local pairing correlation functions (c r |C r+T j) 
are larger than that for the average gap, but still small. Fig. as an example, shows the 
relative error pattern for the local correlation function (c r |C r+T j) for a s-wave pairing state. 
We find that the errors for (c r |C r+T f) are largest (about 5%) at the impurity sites. Clearly 
(c r |C r+T |) is overestimated in the vicinity of impurities and underestimated far away from 
the impurities. Nevertheless it is encouraging that the error made in neglecting off-diagonal 
disorder is smalfl. 

Using the above approximation, we have evaluated the average gaps A for both s- and 
d-wave pairing states in different scattering potentials. For all the cases we have studied, we 
find that A decays almost linearly with x for small x. Fig. |3| shows A r as a function of the 
impurity concentration x for both pairing states in a strong 5-function scattering potential 
Vo = 20 on a 14x14 lattice. All four curve shown in the figure are nearly parallel to each 
other although their values at x=0, i.e. A(0), are quite different. This remarkable behavior 
indicates that the reduction of A by disorder (dA/dx) is determined only by the scattering 
potential to a first approximation and is independent of the pairing symmetries and the 
values of the energy gaps for pure samples. This nearly universal property of the energy 
gap could be of use in future studies of disordered superconductors. For weak ^-function 
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potentials or finite but short ranged potentials, similar results for A have been found. But 
for a weak 5-function potential the slope of the decay of A with x becomes smaller. 

Our calculations for A and p s have been done mostly on lattices with size ranged from 
10x10 to 18x18 sites. In general, A and p s are size dependent. However, for the quantities 
which are physically more interesting, the relative energy gap A(x)/A(0) and the relative 
superfluid density p s (x)/p s (Q), the finite size effect is small. A comparison for A(x)/A(0) 
and p s (x) / p s (0) for the s- and d-wave states with a ^-function potential Vq = 8 on lattices of 
10x10, 14x14, and 18x18 sites is given in Fig. Similar results have been found for other 
impurity potentials. In general, we find that the finite size effect for p s is larger than that 
for A T when A T is small (i.e. large £). Physically this is because that A rT is determined 
only by the local pair correlation function while p s is determined by the current-current 
correlation function at long wavelengths which in turn is determined by the properties of 
the eigenfunction on the whole lattice. 

Fig. [5] shows A(a;)/A(0) and p s (x) / p s (0) as functions of x for the s- and d-wave pairing 
states in three potentials. For a weak ^-function potential, A(x)/A(0) and p s (x)/ p s (0) 
decrease very slowly with x and the difference between the s- and d-wave pairing states is 
also small. With increasing Vq, the difference between these two pairing states increases. 
The disorder has stronger effect on the d-wave state than the s-wave state; in particularly, 
p s falls much faster in the d-wave state than in the s-wave state. This property may be 
useful for distinguishing a d-wave pairing state from a s-wave pairing state in the unitary 
scattering limit from an experimental point of view. But if an off-site scattering potential 
(i.e. V% term) is present, this difference will be eventually reduced. The disorder effect is 
clearly strongly enhanced in a finite range potential than in a ^-function potential. 

Fig. H shows A(Vo, Vi)/A(0, 0) as functions of V for the s- and d-wave pairing states 
with a;=0.02. The qualitative behaviors of A(Vo, Vi)/A(0, 0) are similar for all the cases 
shown in the figure. A(Vo, Vi)/A(0, 0) decreases with Vq for small Vq, but soon becomes 
saturated when Vq surpasses the band width. dA/dx for the s- and d-wave states can be 
very large depending on the value of V±. 
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B. Density of states 



Now we consider the effect of disorder on the density of states of quasiparticles. We 
calculate the density of states using a recursion methodEl. This method addresses the local 
density of states of an infinite lattice. In this method, the density of states p(E) is obtained 
from the imaginary part of the one-particle Green's function G(E). 

p(E) = lim lmG{E + ie). (10) 

Given a starting state |0), the recursion method is defined by recurrence relations 

#|0) = ao|0) + &i|l) (11) 

and 

H\n) = b n \n — 1) + a n \n) + b n+1 \n + 1) (n>0), (12) 

where {|n)} is a set of normalized bases generated automatically from equations ([11]) and 
(|12|). From the a's and b's generated, G(E) can be expressed in a continued fraction form 

G(E) = g . (13) 

(E - oo) 1 72 

(E - a x ) - 



(E-a 2 ) 

In real calculation, this continued fraction is truncated at a certain step and the remainder 
of the continued fraction is replaced by a parameter which is determined such that the error 
is minimized. We truncate the continued fraction at a step when the difference between the 
result obtained at that step and that with 5 more steps is smaller than the error demanded. 
In using the recursion method, the values of A rr and p obtained previously on finite lattices 
will be used. When the impurity concentration is finite, the approximation A rr = A T 
is assumed. For one impurity system, the strict self-consistent solution for A rr on a small 
lattice around the impurity is used, while for the rest part of the lattice A rT are approximated 
by the average value of A rT on the edges of the small lattice. 
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As discussed by Byers et altJ'EI the densities of states in the vicinity of the impurity is 
in principle measurable via the spatial variation of the tunneling conductance around an 
impurity with a scanning-tunneling-microscope study of the surface of a superconductor. 
In continuum space, the local density of states (or the tunneling conductance) around an 
impurity has been calculated by Byers et afl for both the s- and d-wave pairing states and 
by Choi0 for the d-wave pairing state. They find that the density of states for a give energy 
oscillates in space and depends strongly on the anisotropy of the gap parameter. When the 
energy is larger than A, the oscillation is largest in the directions of the gap maxima and 
smallest in the directions of the gap minima. In lattice space, however, we find that their 
results are partly altered. Fig. ^ shows the impurity induced density of states as a function 
of distance from an impurity along two directions for the s- and d-wave pairing states. (Here 
the fully self-consistent isolated impurity result for the gap function on a 21 x 21 lattice 
shown in Fig. [I] is used as an input.) The density of states oscillates in space with an energy 
and direction dependent wavelength in agreement with the results of RefS In the s-wave 
case, the oscillation along the diagonal direction is much larger than that along the x-axis 
direction, in contrast to the isotropic s-wave pairing state in continuum space. Since the 
energy gaps are the same on the Fermi surface for a s-wave states, this difference is purely 
a lattice effect. For the d-wave pairing state, the oscillations along two directions are not so 
different as shown in Refs. ME. The impurity induced density of states decays slightly faster 
along the x-axis direction than along the diagonal direction. 

To compute the density of states with finite doping of impurities, we have used the results 
of the average gap A obtained previously. The fluctuation of A rr in space is ignored in this 
calculation. Fig. ^| shows the density of states for the s- and d-wave pairing states with 
different potentials and x. The main results are summarized as follows: 

(a) For s-wave pairing with weak potential or strong potential with very small x, p has 
no qualitative change with respect to the case without disorder. In particularly, the energy 
gap still exists and is hardly changed by disorder in agreement with Anderson theorem0. 

(b) For the d-wave pairing state with weak scattering potential, the change of p with 
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respect to the case without disorder is small. But p at the Fermi energy Ep becomes finite, in 
consistent with the non-magnetic impurity scattering theory in the Born scattering limit0. 

(c) For the d-wave pairing state with very strong onsite potential, p shows a peak around 
Ep. This results agrees very well with the self-consistent t-matrix theory for the d-wave 
superconductor in the unitary scattering limitH. On the other hand it also shows that 
the approximations made in the self-consistent t-matrix theory, such as ignoring the vertex 
corrections and the energy dependence of the self-energy, are valid for the d-wave state. 

(d) For the d-wave pairing state with a strong on-site potential or both pairing states with 
a finite range potential, p at Ep grows quickly with increasing x and becomes comparable 
with the average density of states at some critical x. For the s-wave pairing state with a 
strong 5-function potential, a finite gap remains when x is smaller than a critical value x c 
within numerical errors. When x>x c , the gap vanishes (but p s and A are non-zero), p at Ef 
is small and increases slowly with increasing x. 

(e) All singularities of p are suppressed by disorder average in these calculations. We 
have not found any evidence for the singular behavior predicted recently by Nersesyan et 
al@ within numerical error. 

III. CONCLUSION 

We have discussed a straightforward numerical technique which allows detailed effects 
of various impurity potentials on superconductors to be investigated with a BCS mean field 
framework. In particular we evaluated the gap parameter, the superfluid density, and the 
density of states for the s-wave and d-wave superconducting states with non-magnetic im- 
purities, as functions of impurity concentrations and scattering potentials. For one impurity 
systems, the local density of states induced by impurity oscillates in space, in agreement 
with known analytic results. For s-wave pairing, the energy gap and the density of states 
are hardly affected by weak disorder (i.e. either weak scatterers or dilute strong scatterers), 
consistent with Anderson theorem. In dilute impurity limit, our results agree well with the 
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self-consistent t-matrix theory in both Born and unitary scattering limits, and in both s- and 
d- pairing states. For the d-wave pairing state, the density of states at Ep becomes finite 
even for weak scattering potential, consistent with the Born scattering theory of the d-wave 
superconductor. For the d-wave pairing state with strong on-site potential, the density of 
states is in good agreement with the self-consistent t-matrix theory for the d-wave super- 
conductor in the unitary limit. For strong scatterers, the energy gap of the s-wave state 
disappears beyond a critical doping level which is sensitive to the range of the impurity 
potential. A finite range potential is shown to have a stronger effect than a short range 
potential in either pairing state. For Zn doped YBaCuO, experiments find that T c varies 
almost linearly with x and drops about 25% for 2% Zn doping^. If we assume the change of 
T c is equivalent to the change of A at zero temperature, we find that dA/dx in a 5-function 
potential is too small to fit quantitatively with experiments even in unitary scattering limit. 
However for a finite range potential, no such difficulty exists. 
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FIGURES 

FIG. 1. Self-consistent gap function for the s-wave superconducting state on a 21x21 lattice. 
Only half of the lattice is shown. The impurity is located at the center of the lattice. 



FIG. 2. Relative errors of the local pairing correlation function ((c r jC r |)^ — (c r jC r |) A )/(c r |C r |) 
(where (c r |C r j) A and (c r j.c r |) ArT are the results obtained with and without approximation, and 
(c r jC r |) is the average of (c r jC r |) A in space) for the s-wave pairing state on a 14x14 lattice with 
3 impurities and a potential Vo = 20 and V% = 0. The impurities are located on the sites where 
the three highest peak emerge. 

FIG. 3. Energy gap A vs x for both pairing states with the coupling constant J=1.5 and J=2.3 
and the scattering potential Vo = 20 and V\ = on a 14x14 lattice. 25 impurity configurations 
are used in disorder average. 

FIG. 4. Normalized energy gap A(x)/A(0) and normalized superfluid density p s (x)/ p s (0) as 
functions of impurity concentration x on the square lattices of 10x10, 14x14, and 18x18 sites with 
Vo = 8 and Vi=0. A(0) is about 0.5 (0.3) for the s- (d-) wave pairing state, which is about 1/20 band 
width. 100 impurity configurations for 10x10 lattice, and more than 25 impurity configurations 
for 14x14 and 18x18 lattices are used in disorder average. 

FIG. 5. Normalized energy gap A(x)/A(0) and normalized superfluid density p s (x)/ p s (0) as 
functions of x for both s-wave and d-wave superconducting states with three different potentials 
on 10x10 lattices. A(0) are the same as for Fig. ||. 100 configurations of impurities are used in 
disorder average. 

FIG. 6. Normalized energy gap A(Vb, Vi)/A(0, 0) as functions of Vo for the s- and d-wave 
pairing states with x=0.02. The values of A(0,0) are the same as A(0) used in Fig. 100 
configurations of impurities are used in disorder average. 
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FIG. 7. The local density of states induced by an impurity as a function of distance r from an im- 
purity along both the x-axis direction (0°) and the diagonal direction (45°) at energy to = 1.5A max . 
V = 2 and V\ = 0. 

FIG. 8. Density of states as functions of energy E for s- and d-wave superconducting states with 
different impurity concentrations x. The parameters A obtained in Fig. || are used in calculations 
here. More than 1000 configurations of impurities are used in disorder average. 
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